#include <iostream>
#include <fstream>
#include <cstdlib>
#include <ctime>
#include <cmath>

#define NUM_NOISE 18
#define SIZE 46

using namespace std;

void apply_lagrange_interp(double[], double[]);
double lagrange_interp(int, int, double, int, double, int, double, int, double);
bool get_4nearest_neighbor(int, double[], int[]);
void print(double[]);
bool make_noise(double[], double[], int[] );
double estimate_relative_error(double[], double[], int[] );
void copy(double[], double[]);

int main(int argc, char const *argv[])
{
	/*
	(int) number_of_signal
	(int) line (int) sample (double) signal[SIZE]
	*/
	ifstream original_signal_file("signal.txt");
	ofstream log_file_avg_re("log_re.txt", ios::app);
 
	srand( time(NULL) );
		
 	int number_of_signal;
	original_signal_file >> number_of_signal;
	
	int line, sample;
	double data;
	double original_signal[SIZE];
	double noised_signal[SIZE];	int noised_position[NUM_NOISE]; 
	double interpolated_signal[SIZE];
	double diff;
	double avg_of_diff = 0.0;
	int number_of_signal_can_noised = number_of_signal;

	for (int i = 0; i < number_of_signal; ++i)
	{	
		diff = 0;
		original_signal_file >> line >> sample;
		cout << "Sample #" << i+1 << "\tline:\t" << line << "\tsample:\t" << sample << endl;
	
		for (int i = 0; i < SIZE; ++i)
		{
			original_signal_file >> data;
			original_signal[i] = data;
		}
		if ( make_noise(original_signal, noised_signal, noised_position) )
		{

			apply_lagrange_interp(noised_signal, interpolated_signal);

			cout << "Original \t";	print(original_signal);
			cout << "Noised   \t";	print(noised_signal);
			cout << "Interpolated   \t";	print(interpolated_signal);
			diff = estimate_relative_error(original_signal, interpolated_signal, noised_position);
			if ( diff < 0 ) cout << "diff < 0; " ; 
			cout << "Relative error\t" << diff << endl << endl;
			avg_of_diff += diff;
		} else {
			cout << "Can not add " << NUM_NOISE << " white noises;"  << endl << endl;
			number_of_signal_can_noised--;
		}
	}
	cout << "Average(RE) of " << number_of_signal_can_noised << " signals\t" << avg_of_diff * 1.0 / number_of_signal_can_noised << endl << endl;
	log_file_avg_re << avg_of_diff*1.0/number_of_signal_can_noised << endl;
	original_signal_file.close();
	
	return 0;
}

bool make_noise(double original_signal[], double noised_signal[], int noised_position[] )
{
	copy(original_signal, noised_signal);

	int position; int near_neigh[4]; bool good_pos2noise;
	for (int i = 0; i < NUM_NOISE; ++i)
	{
		good_pos2noise = false;

		for (int c = 0; c < 46*4 ; c++ ){
			position = rand() % 40 + 2;
			if ( noised_signal[position] != 0 &&  get_4nearest_neighbor(position, noised_signal, near_neigh) )
			{
				good_pos2noise = true;
				noised_position[i] = position;
				noised_signal[position] = 0;
				break;
			}
		}

		if ( !good_pos2noise )
			return false;
	}
	return true;
}

double estimate_relative_error(double original_signal[], double interpolated_signal[], int noised_position[] )
{
	double sum_re = 0; int pos;
	for (int i = 0; i < NUM_NOISE; ++i)
	{
		pos = noised_position[i];
		sum_re += 1.0 * fabs( original_signal[pos] - interpolated_signal[pos]) / fabs( original_signal[pos] ) ;
	}
	return 1.0 * sum_re / NUM_NOISE;
}

void apply_lagrange_interp( double value[], double lagrange_interp_value[])
{
	for (int i = 0; i < SIZE; ++i)
	{
		lagrange_interp_value[i] = value[i];
		if ( value[i] == 0)
		{
			int neighbor[4];
			if ( get_4nearest_neighbor(i, value, neighbor) )
			{
				lagrange_interp_value[i] = lagrange_interp(i, neighbor[0], value[neighbor[0]], neighbor[1], value[neighbor[1]], neighbor[2], value[neighbor[2]], neighbor[3], value[neighbor[3]] );
			}
		}
	}
}

double lagrange_interp(int x, int x1, double y1, int x2, double y2, int x3, double y3, int x4, double y4)
{
	double l1 = 1.0 * (x-x2)/(x1-x2) * (x-x3)/(x1-x3) * (x-x4)/(x1-x4); 
	double l2 = 1.0 * (x-x1)/(x2-x1) * (x-x3)/(x2-x3) * (x-x4)/(x2-x4);
	double l3 = 1.0 * (x-x1)/(x3-x1) * (x-x2)/(x3-x2) * (x-x4)/(x3-x4);
	double l4 = 1.0 * (x-x1)/(x4-x1) * (x-x2)/(x4-x2) * (x-x3)/(x4-x3);
	return y1*l1 + y2*l2 + y3*l3 + y4*l4;
}

bool get_4nearest_neighbor(int x, double value[], int neighbor[])
{
	int n = 0;
	for (int i = 1; i < 4; ++i)
	{
		if ( (x-i) >= 0 && value[x-i] != 0)
		{
			neighbor[n] = x-i;
			n++; if (n >= 4) break;
		}
		if ( (x+i) < SIZE && value[x+i] != 0)
		{
			neighbor[n] = x+i;
			n++; if (n >= 4) break;
		}
	}
	if (n == 4) return true;
	return false;
}

void print(double *value)
{
	for (int i = 0; i < SIZE; ++i)
	{
		cout << value[i] << "\t";
	}
	cout << endl;
}

void copy(double source[], double dest[])
{
	for (int i = 0; i < SIZE; ++i)
	{
		dest[i] = source[i];
	}
}